Linux Cubed Series 7: Sunsite/Linux Cubed Series 7 - Sunsite Vol 1.iso/system/shells/scsh-0.4/scsh-0/scsh-0.4.2/rts/bignum.scm

Text File | 1995-10-13 | 10KB | 337 lines

; Copyright (c) 1993, 1994 Richard Kelsey and Jonathan Rees.  See file COPYING.

; This is file bignum.scm.

; Integer arithmetic

(define-extended-number-type :bignum (:exact-integer)
  (make-bignum sign magnitude)
  bignum?
  (sign bignum-sign)
  (magnitude bignum-magnitude))

(define (integer->bignum m)
  (if (bignum? m)
      m
      (cond ((>= m 0)
         (make-bignum 1 (integer->magnitude m)))
        ((= m least-non-bignum)
         (make-bignum -1 least-non-bignum-magnitude))
        (else
         (make-bignum -1 (integer->magnitude (- 0 m)))))))

;(define (bignum->integer n)             ;For debugging
;  (* (bignum-sign n)
;     (reduce (lambda (d n) (+ d (* n radix)))
;             0
;             (bignum-magnitude n))))

(define (make-integer sign mag)
  (if (> sign 0)
      (if (smaller-magnitude? greatest-non-bignum-magnitude mag)
      (make-bignum sign mag)
      (magnitude->integer mag))
      (if (smaller-magnitude? least-non-bignum-magnitude mag)
      (make-bignum sign mag)
      (if (same-magnitude? mag least-non-bignum-magnitude)
          least-non-bignum
          (- 0 (magnitude->integer mag))))))

; Arithmetic

(define (integer+ m n)
  (let ((m (integer->bignum m))
    (n (integer->bignum n)))
    (let ((m-sign (bignum-sign m))
      (m-mag (bignum-magnitude m))
      (n-sign (bignum-sign n))
      (n-mag (bignum-magnitude n)))
    (if (= m-sign n-sign)
    (make-integer m-sign (add-magnitudes m-mag n-mag))
    (if (smaller-magnitude? m-mag n-mag)
        (make-integer (- 0 m-sign) (subtract-magnitudes n-mag m-mag))
        (make-integer m-sign (subtract-magnitudes m-mag n-mag)))))))

(define (integer- m n)
  (integer+ m (integer-negate n)))

(define (integer-negate m)
  (cond ((bignum? m)
     (make-integer (- 0 (bignum-sign m))
               (bignum-magnitude m)))
    ((= m least-non-bignum)
     (make-bignum 1 least-non-bignum-magnitude))
    (else (- 0 m))))

(define (integer* m n)
  (let ((m (integer->bignum m))
    (n (integer->bignum n)))
    (make-integer (* (bignum-sign m) (bignum-sign n))
          (multiply-magnitudes
           (bignum-magnitude m)
           (bignum-magnitude n)))))

(define (integer-divide m n cont)
  (let ((m (integer->bignum m))
    (n (integer->bignum n)))
    (divide-magnitudes
       (bignum-magnitude m)
       (bignum-magnitude n)
       (lambda (q r)
     (cont (make-integer (* (bignum-sign m) (bignum-sign n)) q)
           (make-integer (bignum-sign m) r))))))

(define (integer-quotient m n)
  (integer-divide m n (lambda (q r) q)))

(define (integer-remainder m n)
  (integer-divide m n (lambda (q r) r)))

(define integer=
  (lambda (m n)
    (let ((m (integer->bignum m))
      (n (integer->bignum n)))
      (and (= (bignum-sign m) (bignum-sign n))
       (same-magnitude? (bignum-magnitude m)
                (bignum-magnitude n))))))

(define integer<
  (lambda (m n)
    (let ((m (integer->bignum m))
      (n (integer->bignum n)))
      (let ((m-sign (bignum-sign m))
        (n-sign (bignum-sign n)))
    (or (< m-sign n-sign)
        (and (= m-sign n-sign)
         (if (< m-sign 0)
             (smaller-magnitude? (bignum-magnitude n)
                     (bignum-magnitude m))
             (smaller-magnitude? (bignum-magnitude m)
                     (bignum-magnitude n)))))))))


; Magnitude (unsigned integer) arithmetic

(define log-radix 14)            ;Cutting it close here...
(define radix (expt 2 log-radix))

(define greatest-non-bignum (+ (expt 2 28) (- (expt 2 28) 1)))
(define least-non-bignum (* (expt 2 28) -2))

(define zero-magnitude '())
(define zero-magnitude? null?)

(define (low-digit m)
  (if (zero-magnitude? m)
      0
      (car m)))

(define (high-digits m)
  (if (zero-magnitude? m)
      m
      (cdr m)))

(define (adjoin-digit d m)
  (if (and (= d 0) (zero-magnitude? m))
      m
      (cons d m)))

(define (integer->magnitude n)
  (if (= n 0)
      zero-magnitude
      (let ((digit (remainder n radix)))
    (adjoin-digit digit
              (integer->magnitude (quotient n radix))))))

(define (magnitude->integer m)
  (if (zero-magnitude? m)
      0
      (+ (low-digit m)
     (* radix (magnitude->integer (high-digits m))))))

(define greatest-non-bignum-magnitude
  (integer->magnitude greatest-non-bignum))

(define least-non-bignum-magnitude
  (adjoin-digit (- 0 (remainder least-non-bignum radix))
        (integer->magnitude
         (- 0 (quotient least-non-bignum radix)))))

; Combine two numbers digitwise using op.

(define (combine-magnitudes m n op)
  (let recur ((m m) (n n) (carry 0))
    (if (and (zero-magnitude? m) (zero-magnitude? n))
    (integer->magnitude carry)
    (let ((result (+ carry (op (low-digit m) (low-digit n)))))
      (let ((q (quotient result radix))
        (r (remainder result radix)))
        (if (< r 0)
        (adjoin-digit (+ r radix)
                  (recur (high-digits m)
                     (high-digits n)
                     (- q 1)))
        (adjoin-digit r
                  (recur (high-digits m)
                     (high-digits n)
                     q))))))))

(define (add-magnitudes m n)
  (combine-magnitudes m n +))

(define (subtract-magnitudes m n)
  (combine-magnitudes m n -))

; Compare

(define same-magnitude? equal?)

(define (smaller-magnitude? m n)
  (let loop ((m m) (n n) (a #f))
    (cond ((zero-magnitude? m)
       (or (not (zero-magnitude? n)) a))
      ((zero-magnitude? n) #f)
      (else
       (loop (high-digits m)
         (high-digits n)
         (or (< (low-digit m) (low-digit n))
             (and (= (low-digit m) (low-digit n)) a)))))))

; Multiply

(define (multiply-magnitudes m n)
  (let recur ((m m) (a zero-magnitude))
    (if (zero-magnitude? m)
    a
    (let ((a (combine-magnitudes a n (lambda (d e)
                       (+ d (* e (low-digit m)))))))
      (adjoin-digit (low-digit a)
            (recur (high-digits m) (high-digits a)))))))

; Divide m/n: find q and r such that m = q*n + r, where 0 <= r < n.
; Oh no... time to get out Knuth...
; The main thing we don't do that Knuth does is to normalize the
; divisor (n) by shifting it left.

(define (divide-magnitudes m n cont)
  (if (zero-magnitude? (high-digits n))
      (divide-by-digit m (low-digit n)
               (lambda (q r)
             (cont q (adjoin-digit r zero-magnitude))))
      (let recur ((m m) (cont cont))
    (if (smaller-magnitude? m n)
        (cont zero-magnitude m)
        (recur
         (high-digits m)
         (lambda (q r)
           ;; 0 <= r < n  and  d < b
           ;; so  b*r + d < b*n.
           (divide-step (adjoin-digit (low-digit m) r)
                n
                (lambda (q1 r)
                  (cont (adjoin-digit q1 q) r)))))))))

; Divide m by n, where  n <= m < b*n, i.e. 1 <= quotient < b.
; E.g.  if  n = 100  then  100 <= m <= 999
;       if  n = 999  then  999 <= m <= 9989

(define (divide-step m n cont)
  (do ((m-high m (high-digits m-high))
       (n-high n (high-digits n-high)))
      ((zero-magnitude? (high-digits (high-digits n-high)))
       ;; Occasionally q^ is one larger than the actual first digit.
       ;; This loop will never iterate more than once.
       (let loop ((q^ (min (guess-quotient-digit m-high n-high)
               (- radix 1))))
     (let ((r (combine-magnitudes m n (lambda (d e)
                        (- d (* e q^))))))
       (if (improper-magnitude? r)
           ;; (begin (write `(addback ,m ,n ,q^ ,r)) (newline) ...)
           (loop (- q^ 1))
           (cont q^ r)))))))

; Compute q such that [m1 m2 m3] = q*[n1 n2] + r with 0 <= r < [n1 n2]
; Can assume b <= [0 n1 n2] <= [m1 m2 m3] <= [n1 n2 b-1]
; Some examples:
;  m / n :  100[1] / 10[02], 099 / 10, 099[1] / 99[0], 999[8] / 99[99]
; Various hacks are possible to improve performance.  In particular, the
; second division can be eliminated if the divisor is normalized.
; See Knuth.
;  [m1 m2] = q0*[n1] + r0
;  [m1 m2 m3] = q0*[n1 n2] + r^
;  r^ = b*r0 + m3 - q0*n2

(define (guess-quotient-digit m n)
  (let ((n1 (low-digit (high-digits n)))
    (n2 (low-digit n))
    (m1 (low-digit (high-digits (high-digits m))))
    (m2 (low-digit (high-digits m)))
    (m3 (low-digit m)))
    (let ((m12 (+ (* m1 radix) m2)))
      (let ((q0 (quotient m12 n1))
        (r0 (remainder m12 n1)))
    (let ((r^ (- (+ (* radix r0) m3) (* q0 n2)))
          (n12 (+ (* n1 radix) n2)))
      (let ((q1 (quotient r^ n12))
        (r1 (remainder r^ n12)))
        (if (> q1 0)
        (begin (display "This should never happen: q1 = ")
               (write q1) (newline)))
        (let ((q (+ q0 q1)))
          (if (< r1 0) (- q 1) q))))))))

(define (improper-magnitude? m)
  (and (not (zero-magnitude? m))
       (or (< (low-digit m) 0)
       (improper-magnitude? (high-digits m)))))

; Special case of division algorithm for single-digit divisor.

(define (divide-by-digit m d cont)
  (if (= d 0)
      (error "integer division by zero" m d)
      (let recur ((m m) (cont cont))
    (if (and (zero-magnitude? (high-digits m))
         (< (low-digit m) d))
        (cont zero-magnitude (low-digit m))
        (recur (high-digits m)
           (lambda (q r)
             (let ((m1 (+ (low-digit m) (* radix r))))
               (cont (adjoin-digit (quotient m1 d) q)
                 (remainder m1 d)))))))))

;(define (divide-test seed)
;  (let ((random (make-random seed)))
;    (let loop ()
;      (let* ((z1 (integer+ (random) (integer* (random) 10000000)))
;             (z2 (integer+ (random) (integer* (random) 10000000)))
;             (n (max z1 z2))
;             (r (min z1 z2))
;             (q (random))
;             (m (integer+ (integer* n q) r)))
;        (if (not (= n r))
;            (integer-divide m n
;                            (lambda (q1 r1)
;                              (if (and (= q q1) (= r r1))
;                                  (begin (display ".")
;                                         (force-output (current-output-port)))
;                                  (error "division error" m n q q1 r r1)))))
;        (loop)))))


; Extend the generic arithmetic operators.

(define-method &integer? ((n :bignum)) #t)
(define-method &exact? ((n :bignum)) #t)

(define-method &+ ((n1 :exact-integer) (n2 :exact-integer)) (integer+ n1 n2))
(define-method &- ((n1 :exact-integer) (n2 :exact-integer)) (integer- n1 n2))
(define-method &* ((n1 :exact-integer) (n2 :exact-integer)) (integer* n1 n2))
(define-method &= ((n1 :exact-integer) (n2 :exact-integer)) (integer= n1 n2))
(define-method &< ((n1 :exact-integer) (n2 :exact-integer)) (integer< n1 n2))

(define-method "ient ((n1 :exact-integer) (n2 :exact-integer))
  (integer-quotient n1 n2))
(define-method &remainder ((n1 :exact-integer) (n2 :exact-integer))
  (integer-remainder n1 n2))